Загрузка необходимых библиотек

library(tidyverse)
library(ggplot2)
library(ggpubr)
library(plotly)
library(rstatix)
library(corrplot)
library(corrr)
library(factoextra)
library(pheatmap)
library(FactoMineR)
library(ggbiplot)
library(tidymodels)
library(embed)

Задание 1.

life_expectancy <- read_rds("life_expectancy_data.rds")

summary(life_expectancy)
##    Country               Year         Gender          Life expectancy
##  Length:195         Min.   :2019   Length:195         Min.   :55.49  
##  Class :character   1st Qu.:2019   Class :character   1st Qu.:70.02  
##  Mode  :character   Median :2019   Mode  :character   Median :77.55  
##                     Mean   :2019                      Mean   :75.52  
##                     3rd Qu.:2019                      3rd Qu.:80.95  
##                     Max.   :2019                      Max.   :88.10  
##   Unemployment    Infant Mortality      GDP                 GNI           
##  Min.   : 0.178   Min.   : 1.40    Min.   :1.884e+08   Min.   :3.754e+08  
##  1st Qu.: 3.735   1st Qu.: 5.35    1st Qu.:1.117e+10   1st Qu.:1.094e+10  
##  Median : 5.960   Median :13.50    Median :3.967e+10   Median :4.009e+10  
##  Mean   : 8.597   Mean   :19.61    Mean   :4.660e+11   Mean   :4.864e+11  
##  3rd Qu.:10.958   3rd Qu.:30.23    3rd Qu.:2.476e+11   3rd Qu.:2.457e+11  
##  Max.   :36.442   Max.   :75.80    Max.   :2.143e+13   Max.   :2.171e+13  
##  Clean fuels and cooking technologies   Per Capita      
##  Min.   :  0.00                       Min.   :   228.2  
##  1st Qu.: 34.50                       1st Qu.:  2165.3  
##  Median : 80.70                       Median :  6624.8  
##  Mean   : 65.98                       Mean   : 16821.0  
##  3rd Qu.:100.00                       3rd Qu.: 19439.7  
##  Max.   :100.00                       Max.   :175813.9  
##  Mortality caused by road traffic injury Tuberculosis Incidence
##  Min.   : 0.00                           Min.   :  0.0         
##  1st Qu.: 8.20                           1st Qu.: 12.0         
##  Median :16.00                           Median : 46.0         
##  Mean   :17.06                           Mean   :103.8         
##  3rd Qu.:24.00                           3rd Qu.:138.5         
##  Max.   :64.60                           Max.   :654.0         
##  DPT Immunization HepB3 Immunization Measles Immunization Hospital beds   
##  Min.   :35.00    Min.   :35.00      Min.   :37.00        Min.   : 0.200  
##  1st Qu.:85.69    1st Qu.:81.31      1st Qu.:84.85        1st Qu.: 1.301  
##  Median :92.00    Median :91.00      Median :92.00        Median : 2.570  
##  Mean   :87.99    Mean   :86.76      Mean   :87.31        Mean   : 2.997  
##  3rd Qu.:97.00    3rd Qu.:96.00      3rd Qu.:96.50        3rd Qu.: 3.773  
##  Max.   :99.00    Max.   :99.00      Max.   :99.00        Max.   :13.710  
##  Basic sanitation services Tuberculosis treatment Urban population
##  Min.   :  8.632           Min.   :  0.00         Min.   : 13.25  
##  1st Qu.: 62.919           1st Qu.: 73.00         1st Qu.: 41.92  
##  Median : 91.144           Median : 82.00         Median : 58.76  
##  Mean   : 77.380           Mean   : 77.57         Mean   : 59.12  
##  3rd Qu.: 98.582           3rd Qu.: 88.00         3rd Qu.: 78.02  
##  Max.   :100.000           Max.   :100.00         Max.   :100.00  
##  Rural population Non-communicable Mortality  Sucide Rate        continent 
##  Min.   : 0.00    Min.   : 4.40              Min.   : 0.300   Africa  :52  
##  1st Qu.:21.98    1st Qu.:11.85              1st Qu.: 2.050   Americas:38  
##  Median :41.24    Median :17.20              Median : 3.500   Asia    :42  
##  Mean   :40.88    Mean   :17.05              Mean   : 4.802   Europe  :48  
##  3rd Qu.:58.08    3rd Qu.:22.10              3rd Qu.: 6.600   Oceania :15  
##  Max.   :86.75    Max.   :43.70              Max.   :30.100

Задание 2.

plot_ly(data = life_expectancy,
        x = ~ `Basic sanitation services`,
        y = ~  `Tuberculosis Incidence`,
        color = ~ continent) %>%
  layout(
    title = 'Отношение охвата базовым санитарным обеспечением и заболеваемости туберкулезом',
    yaxis = list(title = 'Заболеваемость туберкулезом',
                 zeroline = FALSE),
    xaxis = list(title = 'Процент охвата базовым санитарным обеспечением',
                 zeroline = FALSE))

Задание 3.

life_expectancy %>%
  filter(continent %in% c("Africa", "Americas")) %>%
  ggqqplot(x = "Life expectancy", facet.by = "continent")

 

Для сравнения распределений можем использовать параметрический тест.

stat.test <- life_expectancy %>%
  filter(continent %in% c("Africa", "Americas")) %>%
  rstatix::t_test(`Life expectancy` ~ continent) %>% 
  rstatix::add_xy_position(x = "continent")

stat.test
## # A tibble: 1 × 12
##   .y.       group1 group2    n1    n2 statistic    df        p y.position groups
##   <chr>     <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>      <dbl> <name>
## 1 Life exp… Africa Ameri…    52    38     -12.3  85.8 1.31e-20       86.6 <chr> 
## # ℹ 2 more variables: xmin <dbl>, xmax <dbl>
life_expectancy %>%
  filter(continent %in% c("Africa", "Americas")) %>%
  ggboxplot(
  x = "continent", y = "Life expectancy", 
  ylab = "Life expectancy", xlab = "Continent") + 
  labs(subtitle = get_test_label(stat.test, detailed = TRUE)) + 
  stat_pvalue_manual(stat.test, tip.length = 0)

# Задание 4.

life_expectancy_num <- life_expectancy %>%
  select(where(is.numeric), -Year)

life_expectancy_cor <- cor(life_expectancy_num)

corrplot(life_expectancy_cor, method = 'circle',
         tl.cex = 0.7,
         tl.col = "black")

corrplot(life_expectancy_cor, method = 'number',
         order = 'alphabet',
         type = 'lower',
         tl.cex = 0.7,
         number.cex = 0.5,
         tl.col = "black")

Задание 5.

life_expectancy_scaled <- scale(life_expectancy_num) #Стандартизируем значения

life_expectancy_dist <- dist(life_expectancy_scaled) #Создаем матрицу дистанций

life_expectancy_hc <- hclust(d = life_expectancy_dist, 
                        method = "ward.D2") #Высчитываем дендрограмму кластеров

fviz_dend(life_expectancy_hc, 
          cex = 0.1) #Визуализируем

Задание 6.

pheatmap(life_expectancy_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = life_expectancy_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 4,
         cutree_cols = length(colnames(life_expectancy_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

 

Мы построили иерархическую кластеризацию по строкам и столбцам, что позволило выделить группы наблюдений с “наименьшей дистанцией” между ними. Отчетливо выделилась группа очень небольшого размера (вторая сверху) с высокими GDP и GNP, в остальных группах оба показателя имеют низкие значения. Можно выделить группу (первая сверху) с чуть более высокими значениями показателей Tuberculosis incidence, Mortality caused by road traffic injury, Infant Mortality, Non-communicable Mortality. Также можно выделить отдельный кластер - в этой группе (первой сверху) очень низкие значения показателей иммунизации Measles, DPT и HepB3. И эта же группа выделяется по кластеру показателей Clean fuels and cooking technologies, Life expectancy и Basic sanitation services (низкие значения).

Задание 7.

life_expectancy.pca <- prcomp(life_expectancy_num, 
                        scale = T) #Сохраняем стандартизацию
summary(life_expectancy.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion  0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion  0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
##                           PC15    PC16    PC17    PC18      PC19
## Standard deviation     0.34546 0.26941 0.20224 0.06968 1.017e-15
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.000e+00
## Cumulative Proportion  0.99377 0.99759 0.99974 1.00000 1.000e+00
fviz_pca_var(life_expectancy.pca, repel = T, col.var = "contrib")

 

По результатам PCA-анализа первые 2 главные компоненты объясняют 51% вариации данных. 75% варианции объясняют первые 5 главных компонент, а 90% - первые 9 главных компонент. Для графика использованы первые 2 главные компоненты (соответствуют осям графика, в скобках указан % объясненной дисперсии - 39.9% и 11.6% соответственно). По графику можно выделить группы переменных: 1) Measles Immunization, DPT Immunization, HepB3 Immunization; 2) Basic sanitation services, Life expectancy, Clean fuels and cooking technologies; 3) Urban population, Per capita; 4) GDP, GNI; 5) Infant mortality, Tuberculosis incidence; 6) Non-communicable Mortality, Mortality caused by road traffic injury; 7) остальные. Также видно, что некоторые переменные имеют отрицательную корреляцию - стрелки направлены в противоположные стороны и достаточно близко к окружности, т.е. имеют весомый вклад в главные компоненты (Urban population и Rural population). Переменные, которые внесли наибольший вклад с точки зрения вариации в первых 2 главных компонентах, можно посмотреть или на этом графике (наиболее светлые и близкие к окружности, но т.к. переменных достаточно много - тяжело точно определить визуально), или построив отдельный график.

fviz_pca_var(life_expectancy.pca, 
             select.var = list(contrib = 5), 
             col.var = "contrib")

 

Таким образом, 5 переменных, которые внесли наибольший вклад с точки зрения вариации в первых 2 главных компонентах (суммарно объясняют 75% вариации), - Measles Immunization, DPT Immunization, HepB3 Immunization, Life expectancy и Infant mortality.

Задание 8.

p <- ggbiplot(life_expectancy.pca, 
         scale=0,
         groups = as.factor(life_expectancy$continent), 
         ellipse = T,
         alpha = 0.4) +
  theme_minimal()

p

ggplotly(p,
         tooltip = 'color')

 

При использовании ggbiplot график нечитаемый (т.к. накладываются названия показателей). Можно использовать вместо этого fviz_pca_biplot (есть опция repel), но такой график полноценно не переводится в plotly, насколько я поняла, в plotly не реализовано repel. Во всяком случае у меня не получилось:(

fviz_pca_biplot(life_expectancy.pca, repel = T,
                geom = "point",
                habillage = life_expectancy$continent,
                addEllipses = T, ellipse.level = 0.68,
                alpha = 0.4)

 

Или, поскольку переменные накладываются друг на друга, можно показать на графике, например, только 5 колонок с наиболее весомым вкладом в главные компоненты.

p2 <- fviz_pca_biplot(life_expectancy.pca,
                geom = "point",
                select.var = list(contrib = 5),
                addEllipses = T, ellipse.level = 0.68,
                habillage = life_expectancy$continent,
                alpha = 0.5)

library(gginnards)
p2 <- shift_layers(p2, "GeomPoint", shift = +1) #Вывели выше слой с точками, чтобы в plotly  названия континентов отображались при наведении на точки, а не на эллипсы. Наверное, это можно сделать более элегантно, но я не сообразила как :(

ggplotly(p2,
         tooltip = 'color')

Задание 9.

При разделении наблюдений на группы по континентам особенно выделяется группа стран Европы, в которой высокие показатели Life expectancy, Basic sanitation services, Clean fuels and cooking technologies, Hospital beds. А также группа стран Африки с высокими показателями Infant mortality, Tuberculosis incidence. Для остальных континентов большой разброс всех показателей в положительном и отрицательном направлении.

Задание 10.

umap_prep <- recipe(~., data = life_expectancy_num) %>%
  step_normalize(all_predictors()) %>%
  step_umap(all_predictors()) %>%
  prep() %>%
  juice()

umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) + #  # можно добавить раскраску 
  geom_point(aes(color = life_expectancy$continent), 
             alpha = 0.7, size = 2) +
  labs(color = NULL)

 

При использовании алгоритма UMAP мы видим “сгусток” точек для Африки, а также более-менее плотный сгусток для Европы, соответственно можно детально провести анализ внутри этих кластеров. Для других континентов данные разрознены. Особенность анализа с использованием алгоритма UMAP по сравнению с PCA в том, что при UMAP значения сгруппированы в “сгустки”, однако при этом не анализируются колонки, анализ проводится только по строкам.

Задание 11.

# Эксперимент 1
life_expectancy_rand <- life_expectancy_num %>%
  select(-sample(seq_len(ncol(.)), size = 5))

life_expectancy_rand.pca <- prcomp(life_expectancy_rand, 
                        scale = T)

summary(life_expectancy_rand.pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.3954 1.3490 1.18155 1.06451 0.94693 0.89428 0.78720
## Proportion of Variance 0.4098 0.1300 0.09972 0.08094 0.06405 0.05712 0.04426
## Cumulative Proportion  0.4098 0.5398 0.63955 0.72049 0.78454 0.84166 0.88593
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.64330 0.61574 0.56484 0.44105 0.36146 0.34481 0.20232
## Proportion of Variance 0.02956 0.02708 0.02279 0.01389 0.00933 0.00849 0.00292
## Cumulative Proportion  0.91549 0.94257 0.96536 0.97925 0.98858 0.99708 1.00000
ggbiplot(life_expectancy_rand.pca, 
         scale=0,
         groups = as.factor(life_expectancy$continent), 
         ellipse = T,
         alpha = 0.4) +
  theme_minimal()

# Эксперимент 2
life_expectancy_rand <- life_expectancy_num %>%
  select(-sample(seq_len(ncol(.)), size = 5))

life_expectancy_rand.pca <- prcomp(life_expectancy_rand, 
                        scale = T)

summary(life_expectancy_rand.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     2.3527 1.4310 1.3176 1.1460 0.89742 0.81683 0.70708
## Proportion of Variance 0.3954 0.1463 0.1240 0.0938 0.05753 0.04766 0.03571
## Cumulative Proportion  0.3954 0.5416 0.6656 0.7594 0.81697 0.86463 0.90034
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.67611 0.56633 0.51540 0.41746 0.36205 0.20327 0.07114
## Proportion of Variance 0.03265 0.02291 0.01897 0.01245 0.00936 0.00295 0.00036
## Cumulative Proportion  0.93299 0.95590 0.97488 0.98732 0.99669 0.99964 1.00000
ggbiplot(life_expectancy_rand.pca, 
         scale=0,
         groups = as.factor(life_expectancy$continent), 
         ellipse = T,
         alpha = 0.4) +
  theme_minimal()

# Эксперимент 3
life_expectancy_rand <- life_expectancy_num %>%
  select(-sample(seq_len(ncol(.)), size = 5))

life_expectancy_rand.pca <- prcomp(life_expectancy_rand, 
                        scale = T)

summary(life_expectancy_rand.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.5711 1.3874 1.2418 0.98620 0.89694 0.84900 0.66659
## Proportion of Variance 0.4722 0.1375 0.1101 0.06947 0.05746 0.05149 0.03174
## Cumulative Proportion  0.4722 0.6097 0.7198 0.78930 0.84676 0.89825 0.92998
##                            PC8     PC9    PC10    PC11   PC12    PC13      PC14
## Standard deviation     0.62264 0.48326 0.40962 0.33275 0.2749 0.07033 2.786e-16
## Proportion of Variance 0.02769 0.01668 0.01198 0.00791 0.0054 0.00035 0.000e+00
## Cumulative Proportion  0.95768 0.97436 0.98634 0.99425 0.9997 1.00000 1.000e+00
ggbiplot(life_expectancy_rand.pca, 
         scale=0,
         groups = as.factor(life_expectancy$continent), 
         ellipse = T,
         alpha = 0.4) +
  theme_minimal()

 

При удалении 5 столбцов отмечается повышение кумулятивного процента объясненной вариации. Во всех 3 повторениях этот процент разный, однако во всех случаях выше, чем в исходной оценке на полном датасете. Итоговые представления на биплотах также различаются в 3 экспериментах и по сравнению с исходным анализом. При этом в 3 экспериментах биплоты имеют разную форму - эллипсы меняют форму от очень вытянутой до приближенной к окружности, также в некоторых случаях меняется направление эллипсов. Также несколько меняется величина и направление стрелок. Такие различия между 3 экспериментами связаны с тем, что каждый раз удаляются разные колонки, т.е. меняется анализируемый датасет.

Задание 12.

life_expectancy_dummy <- life_expectancy %>%
  mutate(Africa = ifelse(continent == "Africa", 1, 0),
         Oceania = ifelse(continent == "Oceania", 1, 0)) %>%
  select(where(is.numeric), -Year)

life_expectancy_dummy.pca <- prcomp(life_expectancy_dummy, 
                        scale = T)

summary(life_expectancy_dummy.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.8400 1.4851 1.3960 1.28878 1.11384 1.01306 0.9514
## Proportion of Variance 0.3841 0.1050 0.0928 0.07909 0.05908 0.04887 0.0431
## Cumulative Proportion  0.3841 0.4891 0.5819 0.66098 0.72006 0.76893 0.8120
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.92891 0.84328 0.70113 0.67720 0.62108 0.55110 0.46766
## Proportion of Variance 0.04109 0.03386 0.02341 0.02184 0.01837 0.01446 0.01041
## Cumulative Proportion  0.85312 0.88699 0.91040 0.93223 0.95060 0.96506 0.97548
##                           PC15    PC16    PC17    PC18    PC19    PC20
## Standard deviation     0.39263 0.35858 0.34148 0.26517 0.20133 0.06900
## Proportion of Variance 0.00734 0.00612 0.00555 0.00335 0.00193 0.00023
## Cumulative Proportion  0.98282 0.98894 0.99449 0.99784 0.99977 1.00000
##                             PC21
## Standard deviation     9.958e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00
ggbiplot(life_expectancy_dummy.pca, 
         scale=0,
         groups = as.factor(life_expectancy$continent), 
         ellipse = T,
         alpha = 0.4) +
  theme_minimal()

 

Наблюдаем снижение процента объясненной вариации по сравнению с исходным анализом, что объясняется добавлением 2 новых переменных. При этом размеры и направление эллипсов не изменились (или изменились минимально, незаметно при визуальной оценке). Добавление дамми-колонок не совсем корректно, т.к. PCA предназначен в первую очередь для работы с количественными переменными, которые перед этим проходят нормирование (т.е. нормированные значения количественных переменных и бинарные данные не совсем корректно сопоставлять). Можно провести FAMD, который будет сочетать анализ для количественных и категориальных переменных.